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Abstract 

We show that, for any lattice polytope P C R d , the set int(P) n 
VL d (provided it is non-empty) contains a point whose coefficient of 
asymmetry with respect to P is at most 8d ■ (81 + 7) 2 . If, moreover, 
P is a simplex, then this bound can be improved to 8 • (81 + 7) 2d+1 . 

As an application, we deduce new upper bounds on the volume of 
a lattice polytope, given its dimension and the number of sublattice 
points in its interior. 



1 Introduction 

A lattice polytope in M. d is a convex polytope whose vertices are lattice points, 
that is, points in Z d . For an integer I > 1, let Ii(P) = int(P) n VL d be the 
set of interior points of P whose coordinates are integers divisible by I. 

Of course, some points of h(P) can lie 'close' to dP, the boundary of P. 
However, our Theorem ^ shows that, provided Ii(P) ^ 0, there is w £ h{P) 
with 

ca(w,P) < 8d- (8l + 7) 22d+ \ (1) 

where ca(w, P) is the coefficient of asymmetry of P about w: 

. max{A I w + Ay G P} 

ca(w, P) = max — — -. 

|y|=i maxjA w — Ay G P\ 



*This research was carried out during the author's stay at the Technical University, 
Berlin, sponsored by the German Academic Exchange Service (DAAD) and the Rouse 
Ball Travelling Fund of Trinity College, Cambridge. 
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Although the function in the right-hand side of (jlj) is enormous, the main 
point is that it depends only on d and I. 

We prove an inequality of this type for the case of a simplex S first. 
Namely, Theorem Q implies that, for some w £ Ii(S), 

ca(w,5) < 8- (8/ + 7) 2d+1 . (2) 

Here the claim essentially concerns the barycentric coordinates (ao, • • • , ce d ) 
of w inside S because of the easy relation 

ca(w, S) = max = ; — - — 1, (3) 

o<i<d oti m s (w) 

where ms(w) := mino<i<dai is the smallest barycentric coordinate of w £ 
5. Define 

(3(d,l) := inf max{m 5 (w) | w £ Ii(S)}, (4) 

where the infimum is taken over all lattice simplices S with Ii(S) ^ 0. 
(For example, it is easy to see that (3(1,1) = jxjO Thus we have to prove 
a positive lower bound on /3(d,l). The gist the proof is that if we have 
w £ Ii(S) with ms(w) being 'small', then using one approximation lemma 
of Lagarias and Ziegler |J we can 'jump' to another vertex w' £ Ii(S) with 
ms(w') > m,g(w), see Theorem |2[ 

In fact, one result of Lawrence |Q, Theorem 3] implies that f3(d, 1) > 
but does not give any explicit bound, see Section ^ here. 

It would be interesting to know how far our bounds (|l]) and (^) are 
from the best possible values. The best values that we know arise from the 
following family of lattice simplices. 

Define inductively the sequence td,i by i^/ = Z+l and t^+i^i = t\ ; — td,i+l- 
(This sequence appears in Q.) Consider the simplex 

Bd,i ■= conv{t 1} i ei, . . . , t d _ hl e d _i, t d j e d , -e d } c R d , (5) 

where (ei, . . . , e^) is the standard basis. It is not hard to see that Ii(B d j) = 
{11,1(1 — e d )}, cf. H Proposition 2.6]. We have mB dl (H) = Ws d j(Z(l — 
e ci)) = , 2 1 , '■ the vertex 11 = (I,. . . ,1), for example, has barycentric coor- 
dinates 

(I I I t d > I 1 



M,l td-l,l td,l — 1 t di + 1 t d) i — 1 t di + 1 
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One can show that t^i > (/ + 1) 2 +1 for d > 2 by considering Udj = t^i — 1; 
hence 

P(d,l)<- r ^-<l(l + l)- 2d -\ d>2. (6) 
t d,l 1 

Thus (|2|) establishes the correct type of dependence on d and I, although 
the gap between the bounds is huge. Perhaps, Bdi gives the actual value of 
the function 0(d,l) as well as the sharp bound for (|l]). 

To extend Theorem || to a general lattice polytope P C M. d (Theorem |j), 
we try to find a lattice polytope P' C P with few vertices such that Ii(P') 
and a homothetic copy of P' covers P. The latter condition gives an upper 
bound on ca(w, P) in terms of ca(w, P') for w £ int(P'), see Lemma ^, and 
is satisfied if, for example, P' D S, where S C P is a simplex of the maximum 
volume. But to get a non-empty Ii(P') we may have to add as many as d 
extra vertices to S. It is now possible to define our jumps within P' to get 
the required w G h{P')- However, the bound (||) for (i-polytopes that we 
obtain is comparable with that for 2(i-dimensional simplices; we believe that 
we lose here too much but we have not found any better argument. 

Next, we investigate the following problem. Let p(d, k, I) (resp. s(d, k, I)) 
be the maximum volume of a lattice polytope (resp. simplex) P C M. d with 
|/;(P)[ = k. As for any d > 2 there exist lattice simplices with no lat- 
tice points in the interior and of arbitrarily large volume, we restrict our 
consideration to the case k > 1. 

Trivially, p(l,k, I) = s(l,k,l) = (k + 1)/. A result of Scott implies 
that p(2, 1, 1) = s(2, 1, 1) = | and p(2, k, 1) = s(2, k, 1) = 2(Jfe + 1) for k > 2. 
Hensley |||, Theorem 3.6] showed that p(d,k,l) exists (i.e., it is finite) for 
k > 1. The method of Hensley was sharpened by Lagarias and Ziegler j|, 
Theorem 1], who showed that 

p(d,k,l) < kl d {7(kl + l)) d2d+ \ (7) 

and also observed that, for any fixed (d, k, I), there are finitely many (up to 
a GL n (Z)-equivalence) lattice polytopes P cM. d with |J/(P)| = k > 1. 

Lagarias and Ziegler Theorem 2.5] proved the following extension of 
a theorem of Mahler [§: "A convex body K C M d with k = \Ii{K)\ > 1 
satisfies 

vol(K) < (Z(ca(w, if) + l)) d ■ k, (8) 

for any w G h{K).' n 
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Combining (§) with ([[]) ( or more exactly with (p7|)), we obtain that 

p(d, k, I) < (8dl) d ■ (81 + 7) d ' 22d+1 • k. (9) 

A theorem of Blichfeldt []l|] says that |PnZ rf | < n + n! vol(P); combined 
with (|9|) it gives an upper bound on |PnZ rf | in terms of \Ii(P)\ (if the latter 
set is non-empty). 

An upper bound on s(d,k,l) can be obtained by applying (|8|) to (^). 
However, we obtain a better bound in Theorem |6| by exploiting the geometry 
of a simplex, namely we show that 

s(d, k, I) < 2 3d - 2 ■ l d ■ (81 + 7)( d - 1 ) 2d+1 . k/dl. (10) 

The best lower bound on p(d, k, I) and s(d, k, I) that we know (except for 
(d, k, I) = (2, 1, 1)), comes from the consideration of the simplex 

Sd,k,l ■= conv{0, ti t i ei, . . . , t d - lt i e d _i, (k + l){t d ,i - l)e d }, 

which satisfies Ii(Sd,k,i) = {Zl + iZe^ | < i < k — 1}, see [^, Proposition 5.6]. 
This demonstrates that 

s (d,k,l)>vol(S dAl )>^-(l + l) 2d -\ 

see formula (2.13) in ||. The family (Sd,k,i) w as found by Zaks, Perles and 
Wills H and its generalization (the addition of parameter /) — by Lagarias 
and Ziegler ||. 

Again, we have the correct type of dependence of d, k and I but the 
gap between the known bounds is huge. The ultimate aim would be to find 
exact values, which is probably not hopeless because the above contructions, 
believed to be extremal, are rather simple. 

2 Jumping inside a simplex 

We will use the following lemma of Lagarias and Ziegler 0, Lemma 2.1]. 
Lemma 1 For a real A > 1 and integer n > 1, define 

5(n,\) = (7(\ + l)y 2n+1 . (11) 



4 



Then, for all positive real numbers ai,...,a n satisfying 

n 

1 - 5(n, A) < y^Oj < 1, 

i=l 

i/iere exist non-negative integers P\, . . . , P n ,Q such that 

Q = p 1 + ... + p n > 0, (12) 
A-P 

> AQ + 1 ^ or 1 - * - n ' ( 13 ) 

AQ + 1 < ^A)- 1 . | (14) 

The above lemma is the main ingredient in our 'jumps.' Here it is applied 
with A = I / . There is nothing special about the constant | except it 
makes ([ll]) look simpler; any fixed number greater than 1 would do as well. 

Theorem 2 Let I > 1 and let S = convjvo, . . . , v^} C M. n be a lattice 
simplex. If rel-int(S') n VU 1 is non-empty, then it contains a point w with 

m s (w) > 7 := 6(d, § 0/8 = (81 + 7)- 2d+1 /8. (15) 

Proof. We may assume that n = d because we can always find a linear 
transformation preserving the lattice Z n (and so lZ n as well) and mapping 
S into R^CR". 

Let w = Ylt=o a i v i e Ii(S), J2i=o ai = 1) be a vertex maximizing 
ms(w). Suppose that the claim is not true. Assume that ao < • • • < a^; 
then ms(w) = ao < 7. Let j be the index with ay < 87 < ay+i; note that 
j < d — 1 is well-defined as ad > gxj > 87. 

We have X^i=o a i < ^7(7 + 1) w hich, as it is easy to see, does not exceed 
5(d — j, 1 1) for j G [0, d — 1]. Hence, Lemma Q is applicable to the d — j 
numbers oy+i, ■ ■ ■ ,ad and yields integers Pj+i, ■ ■ ■ , -Pd, Q satisfying (12)- 

Consider the vertex 

d 

W ' = (IQ + 1)W - lP M G 

i=j+l 

We have w' = 2~^f=o a i v *' where, for i £ [0, j], := (/Q + l)aj > ao and, 
for i£[j + 1, d], a[ := (IQ + l)ai - IP { > a { /8 > a by (||). As 

n d d 

= (lQ + l)J2 a i~ Yl * P * = 1 ' 

j=0 i=0 i=j+l 
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the lattice point w' lies in the interior of S and contradicts the choice of w. | 

Remark. For n = d the inequality (g) claimed in the introduction follows 
by applying (B) to the vertex w G Ii(S) given by Theorem 



3 I) for small d and / 

Let us try to deduce some estimates of (3(d,l) when d and I are small. We 
have a general upper bound ([]) which in particular says that 

,3(3,1) < 1, 0(3,2) < Jj. (17) 

Here we present some results obtained with the help of computer showing 
that ([l6|) and (|l7| ) are probably sharp. 

How can one get a lower bound on, for example, (3(2, 1)? Our approach 
is the following. 

Given a lattice simplex S = conv{vo, Vi, V2} C M 2 , let w be a lattice 
vertex maximizing ms over I\(S) 7^ 0. Write the barycentric representation 
w = Ya=q a i v i- We nave 

q + «i + a 2 = 1. (18) 

Without loss of generality we may assume that 

< «o < Q!i < a.%. (19) 

Consider the vertex w' = 2w — V3 which is the jump of w corresponding 
to P = (0,0,1). Its barycentric coordinates (ckq, a'i, a^) = (2ao, 2«i, 2«2 — 1) 
satisfy a' x > Qq > ao- By the choice of w, we must have 

2a 2 - 1 < a - (20) 

Similarly, the (0, 1, l)-jump w' = 3w — vi — v 2 satisfies a' = 3«o > 
and, in the view of a' 2 > a' l5 we obtain 

3ai - 1 < a . (21) 
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Not everything goes so smoothly if we consider e.g. the (0, 1, 2)-jump, 
when we can only deduce that 



4qi - 1 < a>o or 4a2 - 2 < ao- (22) 



How small can ao be, given only the constraints (which can 



be realized as a mixed integer program)? Solving this MIP, we obtain that 
ao > to- Knowing this bound, we can enlarge our arsenal of jumps. For 
example, the vertex w' = 12w — vi — 4v2 — 6V2 satisfies a' = 12ao — 1 > «o; 
hence 

12ai - 4 < a or 12a 2 - 6 < a . (23) 



The addition of (|23|) to the system (081)— (22) , improves our lower bound 



to a > Yf. 

In this manner we can repeatedly add new constraints to our MIP as long 
as this improves the lower bound on ao- This was realized as a program in C 
which can be linked with either CPLEX (commercial) or lp_solve (public 
domain) MIP solver. The source code is freely available || and the reader 
is welcome to experiment with it. 

When one runs the program, it seems that the obtained lower bound / 
on ao approaches some limit without attaining it. And, of course, the more 
iterations we perform, the larger coefficients in the added inequalities are and 
the problem becomes more complex. Hence, the question to what extend we 
can trust the output should be considered. The CPLEX has the mechanism 
to set up various tolerances which specify how far CPLEX allows variables to 
violate the bounds and to be still considered feasible during perturbations. 
The manual does not specify the guaranteed accuracy. As the coefficient 
at «o is always 1, we assume that the error does not exceed A, the largest 
sum of absolute values of coefficients in one inequality multiplied by the 
tolerance, which we set to 10 -9 , the smallest value allowed by CPLEX's 
manual. 

In particular, we allow a (Pq, Pi, i^H ump when we can guarantee that 
(1 + P + Pi + P 2 )a -P > a , namely, when (/ - A)(P + P l + P 2 ) > P . 

We ran the program, with CPLEX 6.6, for various d and I; Table || 
records the obtained lower bounds up to 10 -6 . Unfortunately, we had no 
success when d > 4 or when d = 3 and / > 3 or when d = 2 and I > 28: 
the obtained lower bound was still zero when the MIP became too large to 
solve. 
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Table 1: Computed lower bounds on (3(d,l) using CPLEX. 
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4 Extending results to lattice polytopes 



First we have to express analytically the intuitively obvious fact that if two 
polytopes cover each other (up to a small homothety) then their coefficients 
of asymmetry cannot be far apart. 

Lemma 3 Let P' C P be two polytopes such that P can be covered by a 
translate of XP' . Then, for any w £ int(P'), 

ca(w,P) < |A|ca(w,P) + |A| - 1. (24) 

Proof. Assume that |A| > 1, for otherwise P' = P and we are home. Also, 
the case of 1-dimensional polytopes is trivial. 

Let wi,w"2 £ dP be two points with w G [wi,W2] and ca(w, P) = 
|wi — w| : |w — W2|. Let dP' fl [wj, w] = {w^}, i = 1, 2, where [x, y] denotes 
the straight line segment between x and y. Clearly, 

|wi — W2I = (ca(w, P) + 1) |w — W2I > (ca(w, P) + 1) [w — w 2 |. (25) 

As XP' covers {wi,W2}, there are ui,U2 € P' with 

wi - w 2 = |A| (ui - u 2 ). (26) 

We can assume that ui,U2 € dP'. If ui = and U2 = w 2 , then we 
let v = U2. Otherwise let v be the (well-defined) point of intersection of 
L(ui,w) and L(u2,w 2 ), where L(x, y) denotes the line though the points x 
and y. As v £ L(ii2,w 2 ) lies outside of int(P'), we have 

( D /\\l u i _w l |ui - u 2 | - |w - w 2 | 
ca(w, P ) > — 



w — V w 



which implies the required by (25) and (]26|). I 



Remark. Note that the bound in ( |24| ) is sharp, as is demonstrated e.g. by 
C d C |A|Cd and w = cl with < c < |, where C d C R d is the 0/1-cube. 

Now we are ready to prove our result on lattice polytopes. 

Theorem 4 Let I > 1 be an integer and let P C M d be a lattice polytope 
with Ii(P) / 0. Then there is w € I Z (P) with 



Sd 



ca(w ' p)£ ^^ffo- 1 = 8<i - (8,+7f+ '- L (27) 
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Proof. Let S = convjvo, . . . v^} C P be a simplex of the maximum volume; 
we may assume that each Vj is a vertex of P. 

Choose u € Ii(P)- Let ui € int(S') be any vertex and let 112 be the point 
of intersection of the ray {ui + A(u — Ui) A > 0} with the boundary of P. 
The vertex 112 lies in the interior of some face which is spanned by at most 
d vertices of P. Hence u 6 rel-int([ui, 112]) can be represented as a positive 
convex combination of n < 2d + 1 vertices of P including all vertices of S, 
let us say u = J2i=o a i v » W1 th Y^=o a i = ^ an< ^ eacn a 'i > 0. 

Choose a vertex w £ Ii(P') and a representation w = X^=o aiVi W1 ^h 
SILo a * = maximizing mino<i< ra ai- Denote this maximum by m(w) > 
0. The argument of Theorem [2] shows that m(w) > 5{n — 1, = Z)/8 > 
5(2d, f0/8. 

The polytope P' can be represented as a projection of an n-simplex S n 
such that w is the image of v € int(S' n ) with ms n (v) = m(w). Now, it is easy 
to see that a linear mapping cannot increase the coefficient of asymmetry; 
hence 

/\ n 1 — m(w) 

ca(w,P ) < ca(v, S n ) = . 

m(w) 

It is known that P C {—d)S + (d + l)s, where s is the centroid of S, see 
e.g. [||, Theorem 3] . By Lemma |||, we obtain 

ca(w, P) < dca(w, P') + d - 1 < d 1 ~ m(w) +d-l = — ^— - 1, 

m(w) m(w) 



which give the required by (|Tl]). I 

Remark. The bound (^) is much worse than (||); the reason is that we 
may have to approximate 2<i-tuples of numbers in Lemma |l[ Unfortunately, 
we cannot guarantee that P' has much fewer than 2d + 1 vertices, as e.g. 
P = conv{±ei, . . . , ie^-i, (k + l)led} demonstrates. Perhaps, one can show 
that any such example cannot be extremal for our problem and thus improve 

on <m. 



Remark. It should be possible to generalize Theorem ^ and || by proving 
the existence of a number b = b(d,l,m) > such that any lattice polytope 
P contains m distinct points in Ii(P) (provided |ij(P)| > m) with coefficient 
of asymmetry of each being at least b. The idea of the proof is the following. 
If P is a simplex, take distinct wi, . . . , w m E Ii(P) with with largest mp's. 
Now each jump of either does not increase mp(wj) or maps w, into some 
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other w j . We are done if we can show that if mp is very small then there are 
at least m distinct jumps increasing it. The latter claim would be achieved 
by rewriting the proof of Lemma EL so that in the conclusion we have at least 
m suitable (n + l)-tuples of integers. To extend the claim to general lattice 
polytopes, observe that m vertices in i/(P) can be represented each as a 
positive combination of d+ 1 vertices of a max- volume symplex and at most 
md other vertices of P and follow the argument of Theorem |j. We do not 
see any principal difficulties arising here, but it would take too much space 
to write the complete proof, so we restrict ourselves to this little observation 
only. 



5 Volume of lattice simplices 

For simplices we have a better method (than applying (Q)) for bounding 
volume which appears in ||, Theorem 3.4] (see also || Lemma 2.3]). Let us 
reproduce this simple argument here. 

Lemma 5 Let S = convjvo, . . . , v^} be any simplex and let w £ h{S) have 
barycentric coordinates (qo, • • • , ay). Then 

l d 

vol(S) < 11,(5)1. (28) 

d\ x ct\ x a2 x • • • x ad 

Proof. Consider the region 

X = {w + YLx A(vi - v ) | |A| < oii 1 < i < d}. 

It is a centrally symmetric parallelepiped around the vertex w € VL d with 
volume vol(X) = d\ vol(S') nf=i(2 Q ? )- We have to show that the volume of 
X cannot exceed (2l) d \Ii(S)\. If this is not true, then X contains (besides 
w) at least |ij(5)| pairs of vertices w ± u € VL d by Corput's theorem ||] 
and, clearly, at least one vertex of each such pair lies within Ii(S), which is 
a contradiction. | 

Now we can deduce the following result. 
Theorem 6 For any k > 1, the inequality (|l~0|) holds. 
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Proof. Let S C W 1 be lattice simplex with |iz(/S)| = k. By Theorem |2] there 
is w G h{S) with mj(w) > 7 = (8Z + 7)~ 2 Assume ao < • • • < otd- 

Then it is easy to see that n?=i <*i > 7 d_1 (1 - 7<0- Th e claim now follows 
from (H|). I 



6 s(d,k,l) for small d and / 

As we have already mentioned, s(2,k,l) was computed by Scott |^]. The 
simplex S^fc,; shows that s(2, k, I) > 1(1 + l) 2 (k + l)/2. Upper bounds on 



s(2, k, I) can be obtained by applying fl28|) to the lower bounds on (3(2,1) 
from Table |l|. But even if we knew j3(2, 1) exactly, the best upper bound on 
s(2,k,l) that this method would give is l 5 k/2 + 0(l 4 )k, so there would still 
be an uncertainty about s(2, k, I). 

Also, an interesting problem is the determination of s(3, k, 1). The sim- 
plex £3^1 shows that s(3, k, 1) > 6(k + 1). Theorem |6| gives, already for 
such small d, very bad bounds. However, there is a very simple argument, 
following the lines of Section |3| and proving that 

29791 

s(3,M) 14.106- A;. ( 29 ) 

Given a lattice simplex S C M 3 , we can deduce as before that the 
barycentric coordinates (ao, . . . , 03) of a lattice vertex maximizing mg sat- 
isfy 203 — 1 < ao, 3a2 — 1 < «o and either 4«2 — 1 < ao or 4«3 — 2 < ao- 
These inequalities do not guarantee yet that Qo > 0, but they guarantee 
that ai > and, as it is routine to see, that 0110.20.3 > X |j x i|, which 
implies @ by (||). 

Of course, we could write more equations on the a's, but this method 
would not lead to the best possible bound. For example, the simplex 153,1,1 
shows that we cannot guarantee a vertex in Ii(S) with 010203 > \ x | x 
so the best bound we would hope to obtain this way is s(3, k, 1) < 12k only. 



7 Lawrence's Finiteness Theorem 

A result of Lawrence @, Lemma 5] implies that there is 7 = 7(d) > 
such that, for any lattice simplex S = convjvo, . . . , v^}, the set h(S) (if 
non-empty) contains a vertex w with oq > 7, where (oq, . . . , Od) are the 
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barycentric coordinates of w. This directly follows from our Theorem ^ 
but unfortunately I could not find a simple argument giving the converse 
implication. And, in fact, the corresponding extremal functions are different: 
for example, (3(2, 1) < 1/8 while it is claimed in |Q, p. 439] that we can take 
7 (2) = 1/6. 

However, one can deduce from j| that (3(d, 1) > using the follow- 
ing simple modification of Lawrence's proof. For i G N let Ui = {w G 
A d | ca(w, A d ) < i}, where A d = conv{0,ei, . . . ,e d }. Clearly, U^ d+1 Ui = 
int(A^). By 0, Theorem 3] there exists i such that for any w G int(A^) there 
are j G N and u G Z d with with jw + u G Ui. We claim that (3(d,l) > j^. 
Indeed, let S C M rf be any lattice simplex and let v G I\{S). Choose any 
affine function / : R d -> M. d with f(S) = A d and let w = /(v) G int(A d ). 
Given w, choose the corresponding j G N and u G TL d . It is easy to check 
that v' = /~ 1 (jw + u) belongs to h(S) and satisfies ms(v') > j^. 

But as it has already been remarked, Lawrence's argument does not give 
any explicit bound. 
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C source code 



/* 

program for proving lower bounds on beta(d,l), the function defined in 
e-print O.Pikhurko "Lattice Points inside Lattice Polytopes" at arXiv.org 

Copyright (C) 2000 Oleg Pikhurko 

This program is free software; you can redistribute it and/or 
modify it under the terms of the GNU General Public License 
as published by the Free Software Foundation. 

This program is distributed in the hope that it will be useful, 
but WITHOUT ANY WARRANTY; without even the implied warranty of 
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 
GNU General Public License for more details: 

\protect\vrule widthOpt\protect\href {http : //www . gnu . org/copylef t/gpl . htmlHhttp : //www . gnu . org/copylef 



Purpose 

For a lattice simplex $S\subset R~d$ and $x\in S$, let $m(x)$ be the 
value of the smallest of the $d+l$ barycentric coordinates of 
$x$. There is $\beta(d,l)>0$ such that if $S$ contains at least one 
vertex of $lZ~d$ in its interior, then it contains such a vertex $x$ 
with $m(x)\ge \beta(d,l)$. The program computes lower bounds on 
$\beta(d,l)$ for small $d$ and $1$. 



Algorithm 

Main idea: let $x\in S$ be a pair (almost) attaining $\beta(d,l)$. Let 
$ (a_0\le\dots\le a_d)$ be the barycentric coordinates of $x$ . For any 
numbers $p_0 , \dots ,p_d\in N$, there is $i\in[0,d]$ such that 
$(lq+l)a_i-lp_i\le a_0$, where $q=\sum_{i=0}~d p_i$. (For otherwise we 
have a contradiction, the vertex $x'=(q+l)x-\sum_{i=l}~d p_i v_i\in 
lZ~d$ belongs to the interior of $S$ and has larger $m$ , where $v_i$'s 
are the vertices of $S$.) 

That is, the $a_i$'s (for extremal $x$) must satisfy certain 

linear inequalities. If these imply a lower bound on $a_0$ it is also 

a lower bound on $\beta(d,l)$. 

In the algorithm, we define our initial LP to consists of $0\le 
a_0\le\dots\le a_d$ and $\sum_{i=0}"d a_i=l$. Then we repeat the 
following. Given LP, let $(a_i)$ be a solution minimizing $a_0$. Try 
to find $p_i$ with $(lq+l)a_i-lp_i>a_0$ which show that $(a_i)$ cannot 
be the coordinates of extremal $x$, add the corresponding constraints 
to our LP, and repeat. 

Of all $p_i$'s, we choose the smallest in the lex order with minimum 



$q$, which probably speeds convergence. Any better ideas? 

As each time we have or-connected inequalities, we introduce binary 
variables; of course, if $p_i=p_{i-l}$ then there is not need to 
include the $i$th inequality. Also the $0$th equality is never 
included because $(a_i)$ minimizes $a_0$ given LP, hence if 
$(lq+l)a_0-lp_0>a_0$ for this LP, it will be true for any larger LP. 



Usage 

Prepare file "d-l-in. lp(s) " containing the lp program at which you 
wish to start. (If the file does not exist, the program creates the 
default initial lp.) Tip: rename file "d-l-last . lp(s) " as 
"d-l-in. lp(s) " to start at the place, where the program finished last 
time . 

Run the program. Enter d 1 when asked. As each subsequent prompt, 
enter for non-iteractive mode, -1 to finish calculations, -2 to save 
lp to the file "d-l-out . lp(s) " or a positive number which gives the 
number of iterations to do till next prompt . At each termination of 
the program (except by ~C) , the last LP problem is saved to file 
"d-l-last. lp(s) " . 

The program can be linked with either cplex (www.cplex.com) or lp_solve 

(\protect\vrule widthOpt\protect\href {f tp : //ftp . ics . ele . tue . nl/pub/lp_solve/Mf tp : //ftp . ics . ele . tue . n 
unfortunately not very stable. Required libraries: 

cplex, socket, m, nsl (if compiled with -DCPLEX) ; 

liblpk.a, m, 1 (if compiled with the -DLPS) . 



Bugs 

If we compile with -DLPS and read the initial lp from file, then 
change_lp causes segmentation fault. However, everything works fine if 
we compile with -DCPLEX or the initial problem is created by make_lp() . 

*/ 



#include <stdio.h> 
#include <math.h> 
#include <string.h> 

/* library specific declarations; some of our function also depend 

on the library; these are placed at the end of the file */ 
#ifdef CPLEX 

#include <cplex.h> 
#define REAL double 

REAL tolerance=(REAL) le-9; /* min tolerance allowed by CPLEX */ 

REAL delta= (REAL) le-9; /* initial upper bound on possible error on a[0] ; 

for CPLEX it will be recomputed with each iteration */ 



/* extension of all lp files */ 
#define EXT ".lp" 
CPXENVptr env = NULL; 
CPXLPptr lp = NULL; 
int status = 0; 

#elif defined LPS 

#include <lpkit.h> 

/* set different extension to avoid interference */ 
#define EXT " .lps" 
lprec *lp; 

REAL delta=(REAL)le-5; /* this is our guess for lp_solve */ 
#else 

#error Please compile with either -DCPLEX or -DLPS option 
#endif 



#define DMAX 9 /* size allocated for a[] ; 

d must not exceed DMAX */ 
#define QMAX 1000 /* if next_p we get as far as q>=QMAX, 

we assume that no approximation exists */ 
#define IMAX 600 /* the maximal number of iterations */ 

/* most important functions */ 

void get_a(void) , read_a(void) ; /* get/(read from stdin) next solution a[] 
void next_p(void) ; /* compute approximation p [] */ 

void change_lp(void) ; /* add the corresponding constraints to lp */ 
/* auxiliary functions */ 

void print_a(REAL*) , print_p(int*) , print_all (void) ; /* various printing */ 

void die (char*); /* terminate the program */ 

void close_all (void) ; /* close all files, etc */ 

void init_lp(void) ; /* initialization */ 

void add_dl (char*) ; 

void f ind_ jumps (void) ; 

void p_string(char*,int) ; 

void save_lp(char*) ; 

void check_solution(void) ; 

void change_lp_plain(void) ; 

void change_lp_tricky (void) ; 

int d=2,l=l; 

REAL best=(REAL) 0.0; /* the best lower bound we can guarantee so far */ 
REAL e=(REAL) le-6; /* small constant to be used in crucial roundings, 

we search for p [] with (l*q+l)a[i]-l*p[i] > a[0]+e */ 



REAL a[DMAX+l], b[DMAX+l]; /* b [] = (l*q+l) a [] -l*p [] */ 
int p[DMAX+l]; /* a, b & p have d+1 elements each */ 
int q; /* q=p [0] + . . . +p [d] */ 

int jp, jp_ind [DMAX+1] ; /* the number and array of jumps in p[] 
int it; /* the number of iterations */ 

int d_stop; /* gives number of iteration before next prompt; 
if <=0 used for various actions */ 



int mainO 
{ 

int next_stop=0; 

printf ( "Please enter d 1 "); 
scanf 07,d °/.d",&d,&l) ; 

if (d>DMAX I I d<=0 I I 1<=0) 

dieC'Wrong input values for d and 1."); 

init_lp() ; 

for(it=0;it<IMAX;it++) { 
get_a(); /* read_a() ; */ 

check_solution() ; /* can be commented out */ 

if (best<a[0] -delta) /* we can improve our lower bound! */ 
best=a[0] -delta; 

next_p() ; 

print_all() ; 

change_lp() ; 

if (it>=next_stop) { 

printf ("\n\nEnter (non-stop), -1 (end) -2 (save lp)\n") 
printf ("or number (>0) runs till stop "); 
scanf C7.d'\&d_stop) ; 

if (d_stop==-l) 

die ("Exiting. . . ") ; 

if (d_stop==-2) { 

char lp_name [] ="x-x-out" EXT; 
save_lp(lp_name) ; 

} 

if (d_stop==0) 



next_stop=IMAX+l ; 

if (d_stop>0) 

next_stop=it+d_stop ; 

} 

} 

dieC'The maximum number of iterations is exceeded."); 
return 0; 

} 



/* reads a[] from stdout; for testing purposes */ 

void read_a(void) 

{ 

int i ; 
REAL suma; 

printf ( "Enter a\n"); 

for(suma=0,i=0;i<d;i++) { 
scanf ('7,lf &a[i]); 
suma+=a[i] ; 

} 

a[d] =l-suma; 

if (a[0]<0 I I a[d]<0) 

die("a[] is not feasible."); 

} 



/* for given a[] compute a smallest p [] such that min( (lq+1) a[] -l*p [] ) >a[0] . 

return if no reasonably small P exists */ 
void next_p(void) 
{ 

int f ound=0 , i , sump,dif f ; 

/* printf ("next_p: obtained the following vector :\n"); 
print_a(a) ; */ 

q=0; /* we try to find p [] with smallest q, so we interate on q */ 
do { 



q++; 



sump=0 ; 

for(i=0;i<=d;i++) { 

/* the case i==0 receives a special treatment */ 
p[i]=ceil(q*(a[i]-(i?0: delta)) + (a[i] -a[0] )/l - 1 

if (p[i]<0) /* this may happen if a[i] is small */ 

p[i]=0; 
sump+=p [i] ; 

} 

if (sump>=q) { 
f ound=l ; 

if (sump>q) { /* here different p [] are possible; 

we take the lex-smallest p [] */ 

dif f =sump-q; 

for(i=0; i<=d && diff>0; i++) 
if (p[i]<=diff) { 

diff-=p[i] ; 

p[i]=0; 
} else { 

p[i]-=diff ; 

diff=0; 

} 

} 

} 

} while (! found && q< QMAX) ; 



if (found) { 

/* let us run small check; just in case */ 
f or (sump=0, i=0; i<=d; i++) 
sump+=p [i] ; 

if (sump != q) 

die("next_p: sump does not equal q"); 

for(i=0;i<=d;i++) { 

b[i] = (l*q+l)*a[i]-l*p[i] ; 
if( p[i]<0 II ( b[i]<= a[0] && a[i]>0 ) ) 
die("next_p: the sequence p [] is bad"); 

} 



} else 



dieC'Suitable p cannot be found. Success!?"); 

} 



/* creates a string "r_p [0] _ . . . _p [d] _i" */ 

void p_string(char s [] , int i) 

{ 

int j ; 

char temp_name [255] ; 

strcpy(s, "r") ; 
for(j=0;j<=d; j++) { 

sprintf (temp_name , "_°/»d" ,p [j] ) ; 

strcat(s,temp_name) ; 

} 

sprintf (temp_name , "_7,d" , i) ; 
strcat(s,temp_name) ; 

} 



/* modifies our lp */ 
void change_lp(void) 
{ 

f ind_jumps() ; /* compute the jump in p [] */ 

if(!jp) /* a check just in case */ 

dieC'No jumps in p[], which cannot happen at all."); 

change_lp_tricky () ; 

/* change_lp_plain() ; */ 

/* change_lp_plain() produces lp which are longer to solve, but 
whose binary variables are easier to understand (by humans) ; 
we keep the definition of change_lp_plain() for debugging */ 

} 



/* compute the number jp and array jp_ind[] of jumps in p [] */ 

void f ind_ jumps (void) 

{ 

int i,p_prev; 
jp=0; 

p_prev=p[0] ; 



for(i=l;i<=d;i++) 
if (p [i] >p_prev) { 
jp_ind[jp++]=i; 
p_prev=p [i] ; 

} 

> 



/* prints d+1 array of reals */ 

void print_a(REAL s [] ) 

{ 

int i ; 

for(i=0;i<=d;i++) 

printf ( ,,0 /„19.17f ",s[i]); 
printf ("\n") ; 

} 



/* prints d+1 array of integers */ 

void print_p(int s [] ) 

{ 

int i ; 

f or(i=0; i<=d; i++) 

printf ("Xd ",s[i]) ; 
printf ("\n") ; 

} 



/* prints some global information */ 

void print_all() 

{ 



printf ( "\nGlobal variables : \n") ; 

printf ("d=7 d, l=7 d, iteration=7 d, q=7 d, delta=7o . 2e\n" , 
d, 1, it, q, delta) ; 

#ifdef CPLEX 

if (status) { /* get the error message corresponding to status */ 
char errmsg[1024] ; 



CPXgeterrorstring (env, status, errmsg) ; 
printf ("7,s\n", errmsg); 

} 

#endif 

printf ("Best guaranteed lower bound so far = 7 19.17f\n", best); 
printf ("Vectors a[] p [] & b[]:\n"); 
print_a(a) ; 
print_p(p) ; 
print_a(b) ; 

} 



/* stop execution */ 
void die (char s [] ) 
{ 

printf ("\n7„s\n",s) ; 
print_all () ; 

close_all () ; 

exit (1) ; 

} 



/* replaces the 1st and 3d latters in a filename by the values of d and 1 */ 

void add_dl(char s [] ) 
{ 

s[0]='0'+d; 
s[2]='0'+l; 

} 



/* check if new a[] is compatible with p [] ; lp_solve sometimes returns 

an unfeasible solution and check_solution fails - any ideas? */ 
void check_solution(void) 
{ 

int i,found=0; 

f or (i=0; i<=d; i++) 

if ((b[i] = (l*q+l)*a[i]-l*p[i]) <= a[0]+e) 
f ound=l ; 



if (Ifound) 



die("check_solution failed."); 

} 



/* library specific definitions of functions: change_lp_plain() , 
change_lp_tricky() , close_all(), save_lp() */ 



#ifdef CPLEX 



/* Initialization */ 
void init_lp(void) 
{ 

int i ; 

char in_lp [] ="x-x-in" EXT; 
FILE *lp_file; 

add_dl(in_lp) ; 

/* Initialize the CPLEX environment */ 
env=CPXopenCPLEXdevelop (festatus) ; 
if (env == NULL) 

dieC'CPXopenCPLEXdevelop failed. ") ; 

/* sets extra space for columns */ 

status=CPXsetintparam(env,CPX_PARAM_COLGROWTH,DMAX*IMAX) ; 
if (status) 

die ("Cannot set CPX_PARAM_COLGROWTH . ") ; 
/* sets extra space for rows */ 

status=CPXsetintparam(env,CPX_PARAM_ROWGROWTH, (DMAX+1) *IMAX) ; 
if (status) 

die ("Cannot set CPX_PARAM_ROWGROWTH . " ) ; 

/* sets extra space for non-zero elements; 

do we need this, given two settings above? */ 
status=CPXsetintparam(env,CPX_PARAM_NZGRDWTH, (DMAX+1) * (DMAX+1) *IMAX) ; 
if (status) 

die ("Cannot set CPX_PARAM_NZGROWTH . " ) ; 

/* sets integrality tolerance to smallest allowable */ 
status=CPXsetdblparam(env,CPX_PARAM_EPINT, tolerance) ; 
if (status) 



die ("Cannot set CPX_PARAM_EPINT. ") ; 



/* sets optimallity tolerance to smallest allowable */ 
status=CPXsetdblparam(env,CPX_PARAM_EPDPT, tolerance) ; 
if (status) 

die ("Cannot set CPX_PARAM_EPOPT . " ) ; 

/* sets feasibility tolerance to smallest allowable */ 
status=CPXsetdblparam(env , CPX_PARAM_EPRHS .tolerance) ; 
if (status) 

die ("Cannot set CPX_PARAM_EPRHS . " ) ; 

/* create lp environment */ 
lp=CPXcreateprob(env,&status, "plattice") ; 
if (lp == NULL) 

dieC'CPXcreateprob failed."); 

if ((lp_f ile=fopen(in_lp, "r"))==NULL) { 

REAL obj [DMAX+1] ; /* objective function */ 

char colname [DMAX+1] [4] ; /* names of variables, "aO", "al", etc */ 

char* colname_p [DMAX+1] ; 

REAL rhs [DMAX+1] ,rmatval [3*DMAX+1] ; 

int rmatbeg [DMAX+2] , rmatind [3*DMAX+1] ; 

char sense [DMAX+1] ; 

printfC'File %s does not exists . \n" , in_lp) ; 
printf ("Creating the basic initial lp.\n"); 

for(i=0;i<=d;i++) { 
obj [i]= (REAL) 0.0; 

colname_p [i] =strcpy (colname [i] , "aO") ; 
colname [i] [l]+=i; 

} 

obj [0]= (REAL) 1.0; 

/* add real, >=0 (by default) variables a0,al,... */ 
status=CPXnewcols ( env , lp , d+1 , obj , NULL , NULL , NULL , colname_p) ; 
if (status) 

dieC'CPXnewcols failed."); 



for(i=0;i<d;i++) { /* inequality -a[i] + a[i+l] >= */ 
rhs [i]= (REAL) 0.0; 
sense [i]='G' ; 
rmatbeg [i] =2*1 ; 
rmatind [2*i] =i ; 
rmatval [2*i] = (REAL) -1 . ; 
rmatind [2*i+l] =i+l ; 
rmatval [2*i+l] =(REAL) 1.0; 

} 



/* equation aO+al+...=l */ 
rhs [d]= (REAL) 1.0; 
sense [d]='E' ; 
rmatbeg[d]=2*d; 
for(i=0;i<=d;i++) { 

rmatind [2*d+i] =i ; 

rmatval [2*d+i] =(REAL) 1.0; 

} 

/* add these basic constraints */ 
status=CPXaddrows (env , lp , , d+1 , 3*d+l , rhs , sense , 

rmatbeg , rmat ind , rmatval , NULL , NULL) ; 

if (status) 

die("CPXaddrows failed."); 

} else { 

f close (lp_f ile) ; 

status=CPXreadcopyprob(env,lp,in_lp, NULL) ; 
if (status) 

dieC'CPXreadcopyprob failed. ") ; 



/* compute a[] as the next solution; */ 

void get_a(void) 

{ 



/* solve the lp */ 

status=CPXmipopt(env, lp) ; 
if (! status) { 

status = CPXgetmipx(env, lp, a, 0, d) ; 
if (status) 

dieC'CPXgetmipx failed\n"); 

> else { /* CPXmipopt failed; maybe there are no bin variables yet */ 

/* printf ("°/.d\n" .status) ; */ 

status=CPXprimopt(env,lp) ; 
if (status) 

dieC'Both CPXmipopt and CPXprimopt failed."); 
status=CPXgetx(env,lp,a,0,d) ; 
if (status) 

dieC'CPXgetx failed. ") ; 



> 

> 



/* add constraints corresponding to p [] to the lp; 

a straightforward version */ 
void change_lp_plain(void) 
{ 

int j , i ; 

int cur_numcols ; /* the number of columns (variables) in lp */ 
/* variables for CPXaddcols */ 

REAL ub [DMAX] ; 

/* variables for CPXaddrows; 

we may need as much as d+1 equations */ 
REAL rhs [DMAX+1] ; 
int rmatbeg [DMAX+2] ; 
int rmatind[4*DMAX] ; 

char sense [DMAX] ; /* types of constraints */ 
char row.name [DMAX+1] [255] ; 

char* row_pointer [DMAX] ; /* symbolic names of rows */ 

REAL rmatval [4*DMAX] ; /* our equations have <=4*d non-zero coefficients */ 
char ctype [DMAX] ; /* types of variables */ 

/* first we add jp binary variables z [0] , . . . ,z [jp-1] */ 
/* for j=l we add an unecessary binary variable, 
but it is fixed to 1 anyway and brings no harm; 
we must also specify ub[j]=l because otherwise 
CPXmipopt returns 3002. */ 
for(j=0;j<jp;j++) { 
ctype [j] = 'B' ; 
lib [j] = (REAL) 1.0; 

} 

status=CPXnewcols (env , lp , jp , NULL , NULL , ub , ctype , NULL) ; 
if (status) 

dieC'CPXnewcols failed."); 

/* get the number of columns (variables) in lp */ 
cur_numcols = CPXgetnumcols (env, lp) ; 
if (cur_numcols<=0) 

die ("CPXgetnumcols failed."); 

/* now out added variables have numbers 
from cur_numcols-jp to cur_numcols-l */ 

/* now we add jp constraints of the form 

-a[0] + (l*q+l)*a[i] + (l*q+l) z[j]<= l*p [i] +(l*q+l) , where i=jp_ind[j] */ 



if (delta<tolerance*(3+2*l*q) ) /* recompute the max abs mistake */ 
delta=tolerance*(3+2*l*q) ; 

for(j=0;j<jp;j++) { 

i=jp_ind[j] ; 

rhs [j] =(REAL) l*p[i]+(l*q+l) ; 

rmatbeg [j] =3* j ; /* three non-zero coeff. per row */ 

rmatval[3*j]=(REAL) -1; 

rmatind [3* j]=0; /* x[0] */ 

rmatval[3*j+l]=(REAL) l*q+l; 

rmatind [3*j+l]=i; /* x[i] */ 

rmatval[3*j+2]=(REAL) l*q+l; 

rmatind [3*j+2] =cur_numcols-jp+j ; /* z[j] */ 

sense [j]='L' ; 

p_string(row_name [j] ,i) ; 

row_pointer [j] =row_name [j] ; 

} 

/* plus the constraint z [0] + . . . +z [jp-1] >=1 ; */ 

if (delta<tolerance*jp) /* recompute the max abs mistake */ 
delta=tolerance* jp ; 

rhs [jp]= (REAL) 1.0; 
sense[jp]='G' ; 
rmatbeg [jp] =3*jp; 
rmatbeg [jp+1] =4*jp; 
p_string(row_name [jp] , 0) ; 
row_pointer [j] =row_name [j] ; 

for(j=0;j<jp;j++) { 

rmatval[3*jp+j]=(REAL) 1; 
rmatind [3*jp+j] =cur_numcols-jp+j ; 

} 

status=CPXaddrows (env , lp , , jp+1 , 4* jp , rhs , sense , 

rmatbeg , rmat ind , rmat val , NULL , r ow_po int er ) ; 

if (status) 

dieO'CPXaddrows failed."); 

} 



/* add constraints corresponding to p[] to the lp; uses the minimum 

nnumber of binary variables; a bit tricky */ 
void change_lp_tricky (void) 
{ 



/* variables for CPXaddrows; 

we may need as much as d equations */ 
REAL rhs [DMAX] ; 
int rmatbeg[DMAX+l] ; 

int rmatind[DMAX*(DMAX+2)] ; /* non-zero elements in added equations */ 
char sense [DMAX] ; /* types of constraints */ 
char row.name [DMAX] [255] ; 

char* row_pointer [DMAX] ; /* symbolic names of rows */ 
REAL rmatval [DMAX* (DMAX+2)] ; 

int j» jp_log=0, j2=l; 

int cur_numcols; /* the number of columns (variables) in lp */ 

/* we add constraints saying that there is j<jp such that 
(l*q+l)a[i]- a[0]<= l*p[i], where i=jp_ind[j]; 
to do so we introduce jp_log binary variables, where 
jp_log is the smallest number with $2~jp_log>= jp. 

For each of 2~jp_log possible combinations of binaries all 
inequalities are vacuous except one, which is precisely one of 
our desired inequalities. We do not explain how such a system 
is constructed but encourage the reader to check the algorithm 
and to run the program and look at the resulting lp */ 

while(j2 < jp) { 
j2*=2; 
jp_log++; 

} 

if (delta<tolerance*(l*q+l)*(jp_log+l)) /* recompute the max mistake */ 
delta=tolerance* (l*q+l) * ( jp_log+l) ; 

if (jp_log>0) { 

/* variables for CPXaddcols */ 
REAL ub [DMAX] ; 

char ctype [DMAX] ; /* types of variables */ 

/* first we add jp_log binary variables 

we must also specify ub[j]=l because otherwise 
CPXmipopt returns 3002. */ 

for(j=0; j<jp_log;j++) { 
ctype [j] = 'B> ; 
ub[j]=(REAL)1.0; 

} 

st atus=CPXnewcols ( env , lp , j p_log , NULL , NULL , ub , ctype , NULL) ; 
if (status) 

dieC'CPXnewcols failed."); 



> 



/* now we add constraints, the case jp_log==0 inclusive */ 

/* get the number of columns (variables) in lp */ 
cur_numcols = CPXgetnumcols (env, lp) ; 
if (cur_numcols<=0) 

die ("CPXgetnumcols failed."); 
/* our added variables have numbers 

from cur_numcols-jp_log to cur_numcols-l */ 

for(j=0;j<jp;j++) { 

int j j , j j 2 , i ; 

i=jp_ind[j] ; 

sense [j]='L' ; 

rmatbeg[j]=j*(2+jp_log) ; /* 2+jp_log non-zero entries per row */ 

rmatind[j*(2+jp_log)]=0; /* coefficient at aO */ 
rmatval [j * (2+jp_log) ] =(REAL) -1.0; 

rmatind[j*(2+jp_log)+l]=i; /* coefficient at a[i] */ 
rmatval [j*(2+jp_log)+l]= (REAL) (l*q+l) ; 

p_string(row_name [j] ,i) ; 
row_pointer [j] =row_name [j] ; 

rhs [j] =(REAL) l*p[i] ; 

for(jj=0,jj2=l; jj<jp_log; jj++,jj2*=2) { 



/* coefficient at jj'th binary variable */ 
rmatind [j * (2+jp_log) +2+j j] =cur_numcols-jp_log+j j ; 



if( ( (unsigned int) 1 « jj) & j ) { 
rhs [j]+=(REAL)l*q+l; 

rmatval [j*(2+jp_log)+2+j j] =(REAL)+l*q+l ; 
} else { 

if (j+jj2>=jp) 

rmatval [ j * (2+ jp_log) +2+ j j ] = (REAL) 0.0; 
else 

rmatval [ j * (2+ jp_log) +2+ j j ] = (REAL) -l*q-l ; 

} 

} 



status=CPXaddrows (env , lp , , jp , jp* (2+ jp_log) , rhs , sense , 



rmatbeg , rmat ind , rmat val , NULL , r ow_po int er ) ; 

if (status) 

dieO'CPXaddrows failed."); 

} 



/* close log-files, cplex and write the current program to file */ 

void close_all(void) 
{ 

char out_lp[]="x-x-last" EXT; 
save_lp(out_lp) ; 

CPXwriteprob(env,lp,out_lp, NULL) ; 
CPXf reeprob(env, &lp) ; 
CPXcloseCPLEX(feenv) ; 

} 



/* saves the lp to file; we cannot here call die() 
because save_lp may be called from it */ 

void save_lp(char s [] ) 
{ 

add_dl(s) ; 

printf ( "Writing current LP to file °/ s ... ",s); 
status=CPXwriteprob(env,lp,s, NULL) ; 

if (status) 

printf ("failed, status=7,d. \n" , status); 
else 

printf ("ok. \n") ; 

} 



#elif defined LPS 



/* read/create the initial lp */ 

void init_lp(void) 

{ 



nstring col_name; 

char in_lp[]="x-x-in" EXT; 



FILE *lp_file; 



add_dl(in_lp) ; 

if ((lp_f ile=f open(in_lp, "r") )==NULL) { 
int i , j ; 

REAL row [DMAX+2] ; 

printfC'File °/.s does not exists . \n" , in_lp) ; 
printf ("Creating the basic initial lp.\n"); 

lp=make_lp(0,d+l) ; 
if (lp==NULL) 

die("make_lp failed."); 

for(i=0;i<=d;i++) { /* give variables names a0,a2,... */ 
strcpy(col_name, "aO") ; 
col_name [1] +=i ; 

set_col_name (lp , i+1 , col_name) ; 

} 

/* constraint -a [i-1] +a [i] >=0 */ 

/* add_constraint starts indexing with 1 and sometimes 
the Oth entry is also used for other purposes; so we 
have to modify all indexing accordingly */ 

for(i=l;i<=d;i++) { 
for(j=0; j<=d+l; j++) 

row [j]= (REAL) 0.0; 
row [i]= (REAL) -1.0; 
row [i+1] = (REAL) 1.0; 

add_constraint(lp,row,GE, (REAL) 0.0) ; 

} 

/* a[0] + . . .+a[d]=l */ 
for(j=l;j<=d+l;j++) 
row [j]= (REAL) 1.0; 
add_constraint(lp,row,EQ, (REAL) 1.0); 

/* set objective function */ 

row [1]= (REAL) 1.0; 
for(j=2;j<=d+l;j++) 
row [j]= (REAL) 0.0; 

set_obj_fn(lp,row) ; 

print_lp(lp) ; 

} else { 

/* NB: keep lp_file open */ 



lp=read_lp_f ile (lp_file , FALSE, "plattice") ; 
f close(lp_f ile) ; 

if (lp==NULL) 

die("read_lp_f ile failed."); 

} 

} 



/* compute a[] as the first d+1 variables of the solution; */ 

void get_a(void) 

{ 

int i ; 

if (solve(lp) ) 

dieC'solve failed."); 

f or(i=0; i<=d; i++) 

a [i] =lp->best_solution [lp->rows+l+i] ; 

} 



/* add constraints corresponding to p [] to the lp; 

a straightforward version */ 
void change_lp_plain(void) 
{ 

nstring row_name; /* ! lp_solve allows at most 24 characters *. 
int j,i,jj; 

REAL col[(DMAX+l)*(IMAX+l)] ; 

/* (DMAX+1) * (IMAX+1) is max possible number of rows/cols in lp 

f or(j=0; j<=lp->rows; j++) 
col [j]= (REAL) 0.0; 

/* first we add jp binary variables z [0] , . . . ,z [jp-1] */ 
/* for j=l we add an unecessary binary variable, 

but it is fixed to 1 anyway and brings no harm; */ 

for(j=0;j<jp;j++) { 
add_column(lp, col) ; 

set_int(lp,lp->columns,TRUE) ; /* the variable is integer */ 
set_upbo (lp , lp->columns , (REAL) 1.0); 
set_lowbo (lp , lp->columns , (REAL) 0.0); 

} 

/* now we add jp constraints of the form 



-a[0]+(l*q+l)a[i]+(l*q+l) z[j]<= l*p [i] +(l*q+l) , where i=jp_ind[j] */ 
for(j=0; j<jp; j++) { 
i=jp_ind[j] ; 

f or ( j j =0 ; j j <=lp->columns ; j j ++) 

col[jj]=(REAL) 0.0; 
col [1]= (REAL) -1.0; /* at a[0] */ 
col [i+l]= (REAL) l*q+l; /* at a[i] */ 
col[lp->columns-jp+j+l]= (REAL) l*q+l; /* at z[j] */ 
add_constraint(lp,col,LE, (REAL) l*p[i]+(l*q+l)) ; 
p_string(row_name,i) ; 
set_row_name (lp , lp->rows , row.name) ; 

} 

/* plus the constraint z [0] + . . . +z [jp-1] >=1 */ 
f or ( j j =0 ; j j <=lp->columns ; j j ++) 

col[jj]=(REAL) 0.0; 
for(j=l; j<=jp; j++) 

col [lp->columns-jp+j] =(REAL) 1.0; 
add_constraint(lp,col,GE, (REAL) 1.0) ; 
p_string(row_name,0) ; 
set_row_name (lp , lp->rows , row_name) ; 

} 



/* add constraints corresponding to p [] to the lp; uses the minimum 

nnumber of binary variables; a bit tricky */ 
void change_lp_tricky (void) 
{ 

int j. jP-l°g=0, j2=l; 

nstring row.name ; /* ! lp_solve allows at most 24 characters */ 
REAL col[(DMAX+l)*(IMAX+l)] ; 

/* (DMAX+1) * (IMAX+1) is max possible number of rows/cols in lp */ 

f or(j=0; j<=lp->rows; j++) 
col [j]= (REAL) 0.0; 

/* we add constraints saying that there is j<jp such that 
(l*q+l)a[i]- a[0]<= l*p[i], where i=jp_ind[j]; 
to do so we introduce jp_log binary variables, where 
jp_log is the smallest number with $2~jp_log>= jp. 

For each of 2~jp_log possible combinations of binaries all 
inequalities are vacuous except one, which is precisely one of 
our desired inequalities. We do not explain how such a system 
is constructed but encourage the reader to check the algorithm 
and to run the program and look at the resulting lp */ 



while(j2 < jp) { 
j2*=2; 



jp_log++; 



/* first we add jp_log binary variables */ 
for(j=0; j<jp_log; j++) { 
add_column(lp, col) ; 

set_int(lp,lp->columns,TRUE) ; /* the variable is integer */ 
set_upbo (lp , lp->columns , (REAL) 1.0); 
set_lowbo (lp , lp->columns , (REAL) 0.0); 



/* our added variables have numbers 

from lp->columns-jp_log to lp->columns-l */ 

/* add constraints now */ 

for(j=0;j<jp;j++) { 

int jj > j j2,i, j j_beg; 
REAL rhs; 

f or ( j j =0 ; j j <=lp->columns ; j j ++) 
col [jj]= (REAL) 0.0; 

i=jp_ind[j] ; 

col [1]= (REAL) -1.0; /* coefficient at aO */ 

col [i+l]= (REAL) l*q+l; /* coefficient at a[i] */ 

rhs=(REAL) l*p[i]; /* rhs as it were without binary variables */ 

for(jj=0,jj2=l; jj<jp_log; jj++,jj2*=2) { 

/* coefficient at jj'th binary variable */ 

j j_beg=lp->columns-jp_log+j j+1 ; 

if( ( (unsigned int) 1 « jj) & j ) { 

rhs+=(REAL)l*q+l; 

col[jj_beg]=(REAL) l*q+l; 
} else { 

if (j+jj2>=jp) 

col [jj_beg]= (REAL) 0.0; 

else 

col [jj_beg]= (REAL) -l*q-l; 

} 

} 

add_constraint (lp, col, LE, rhs); 

p_string(row_name, j) ; 

set_row_name (lp , lp->rows , row_name) ; 



/* write out d-l-last . lp(s) */ 

void close_all(void) 

{ 

char out_lp[]="x-x-last" EXT; 
save_lp(out_lp) ; 

} 



/* saves the lp to file; we cannot here call die() 
because save_lp may be called from it */ 

void save_lp(char s [] ) 



FILE *lp_file; 
add_dl(s) ; 

printf ("Writing current LP to file 7,s ... ",s); 
if ( (lp_f ile=f open(s , "w") )==NULL) 

printf ("failed. \n") ; 
else { 

write_LP(lp,lp_f ile) ; /* purify complains here about "uninitialized memory 

read" , though the code seems clean to me */ 

if (f close (lp_f ile)) 

printf ("failed. \n") ; 
else 

printf ("ok\n") ; 

} 



#endif 



